Build manifold for Neuronal populations¶
- Carlos Talavera-López Ph.D | Wellcome Sanger Institute
- v200217
Import required modules¶
In [2]:
import scgen
import anndata
import numpy as np
import scanpy as sc
import pandas as pd
import tensorflow as tf
In [3]:
sc.settings.verbosity = 3
sc.logging.print_version_and_date()
sc.settings.set_figure_params(dpi = 260, color_map = 'RdPu', dpi_save = 260, vector_friendly = True, format = 'svg')
Fixed function to correct batches.¶
- This was modified for
scGennot to jumble theadataindex and to keepadata.raw. - Once
scGenupdates this in the new version, we can skip it.
In [4]:
import os
from random import shuffle
import anndata
import numpy as np
import scanpy as sc
from matplotlib import pyplot as plt
from scipy import sparse
from sklearn import preprocessing
import pandas as pd
import scgen
def batch_removal_ct5(network, adata, batch_key="batch", cell_label_key="cell_type"):
"""
Removes batch effect of adata
# Parameters
network: `scgen VAE`
Variational Auto-encoder class object after training the network.
adata: `~anndata.AnnData`
Annotated data matrix. adata must have `batch_key` and `cell_label_key` which you pass to the function
in its obs.
# Returns
corrected: `~anndata.AnnData`
Annotated matrix of corrected data consisting of all cell types whether they have batch effect or not.
# Example
```python
import scgen
import anndata
train = anndata.read("data/pancreas.h5ad")
train.obs["cell_type"] = train.obs["celltype"].tolist()
network = scgen.VAEArith(x_dimension=train.shape[1], model_path="./models/batch")
network.train(train_data=train, n_epochs=20)
corrected_adata = scgen.batch_removal(network, train)
```
"""
if sparse.issparse(adata.X):
latent_all = network.to_latent(adata.X.A)
else:
latent_all = network.to_latent(adata.X)
adata_latent = anndata.AnnData(latent_all)
adata_latent.obs = adata.obs.copy(deep=True)
unique_cell_types = np.unique(adata_latent.obs[cell_label_key])
shared_ct = []
not_shared_ct = []
for cell_type in unique_cell_types:
temp_cell = adata_latent[adata_latent.obs[cell_label_key] == cell_type]
if len(np.unique(temp_cell.obs[batch_key])) < 2:
cell_type_ann = adata_latent[adata_latent.obs[cell_label_key] == cell_type]
not_shared_ct.append(cell_type_ann)
continue
temp_cell = adata_latent[adata_latent.obs[cell_label_key] == cell_type]
batch_list = {}
batch_ind = {}
max_batch = 0
max_batch_ind = ""
batches = np.unique(temp_cell.obs[batch_key])
for i in batches:
temp = temp_cell[temp_cell.obs[batch_key] == i]
temp_ind = temp_cell.obs[batch_key] == i
if max_batch < len(temp):
max_batch = len(temp)
max_batch_ind = i
batch_list[i] = temp
batch_ind[i] = temp_ind
max_batch_ann = batch_list[max_batch_ind]
for study in batch_list:
delta = np.average(max_batch_ann.X, axis=0) - np.average(batch_list[study].X, axis=0)
batch_list[study].X = delta + batch_list[study].X
temp_cell[batch_ind[study]].X = batch_list[study].X
shared_ct.append(temp_cell)
all_shared_ann = anndata.AnnData.concatenate(*shared_ct, batch_key="concat_batch", index_unique=None)
if "concat_batch" in all_shared_ann.obs.columns:
del all_shared_ann.obs["concat_batch"]
if len(not_shared_ct) < 1:
corrected = anndata.AnnData(network.reconstruct(all_shared_ann.X, use_data=True))
corrected.obs = all_shared_ann.obs.copy(deep=True)
corrected.var_names = adata.var_names.tolist()
corrected = corrected[adata.obs_names]
if adata.raw is not None:
adata_raw = anndata.AnnData(X=adata.raw.X, var=adata.raw.var)
adata_raw.obs_names = adata.obs_names
corrected.raw = adata_raw
return corrected
else:
all_not_shared_ann = anndata.AnnData.concatenate(*not_shared_ct, batch_key="concat_batch", index_unique=None)
all_corrected_data = anndata.AnnData.concatenate(all_shared_ann, all_not_shared_ann, batch_key="concat_batch", index_unique=None)
if "concat_batch" in all_shared_ann.obs.columns:
del all_corrected_data.obs["concat_batch"]
corrected = anndata.AnnData(network.reconstruct(all_corrected_data.X, use_data=True), )
corrected.obs = pd.concat([all_shared_ann.obs, all_not_shared_ann.obs])
corrected.var_names = adata.var_names.tolist()
corrected = corrected[adata.obs_names]
if adata.raw is not None:
adata_raw = anndata.AnnData(X=adata.raw.X, var=adata.raw.var)
adata_raw.obs_names = adata.obs_names
corrected.raw = adata_raw
return corrected
Read in sources¶
In [5]:
neuronal_cells = sc.read_h5ad('/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/1-source_annotation/heart_neuronal_cells_annotated_ctl200217.h5ad')
neuronal_cd45 = sc.read_h5ad('/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/1-source_annotation/heart_neuronal_cd45_annotated_ctl200217.h5ad')
neuronal_sanger = sc.read_h5ad('/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/1-source_annotation/heart_sanger_neuronal_nuclei_annotated_ctl200217.h5ad')
neuronal_harvard = sc.read_h5ad('/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/1-source_annotation/heart_harvard_neuronal_nuclei_annotated_ctl200217.h5ad')
Change labels for each source¶
In [6]:
neuronal_cells.obs['cell_states'].cat.categories
Out[6]:
In [7]:
trans_from=[['CD9+N0'],['Str1']]
trans_to = ['CD9_Neuronal', 'Stromal-like']
neuronal_cells.obs['cell_type'] = [str(i) for i in neuronal_cells.obs['cell_states']]
for leiden,celltype in zip(trans_from, trans_to):
for leiden_from in leiden:
neuronal_cells.obs['cell_type'][neuronal_cells.obs['cell_type'] == leiden_from] = celltype
In [8]:
neuronal_cd45.obs['cell_states'].cat.categories
Out[8]:
In [9]:
trans_from=[['CD9+N0'],['Str1']]
trans_to = ['CD9_Neuronal', 'Stromal-like']
neuronal_cd45.obs['cell_type'] = [str(i) for i in neuronal_cd45.obs['cell_states']]
for leiden,celltype in zip(trans_from, trans_to):
for leiden_from in leiden:
neuronal_cd45.obs['cell_type'][neuronal_cd45.obs['cell_type'] == leiden_from] = celltype
In [10]:
neuronal_sanger.obs['cell_states'].cat.categories
Out[10]:
In [13]:
trans_from=[['NRXN1_N0'],['CD9+N1'],['Str2'],['doublets_N3'],['N4']]
trans_to = ['NRXN_Neuronal', 'CD9_Neuronal', 'Stromal-like', 'doublets', 'basic']
neuronal_sanger.obs['cell_type'] = [str(i) for i in neuronal_sanger.obs['cell_states']]
for leiden,celltype in zip(trans_from, trans_to):
for leiden_from in leiden:
neuronal_sanger.obs['cell_type'][neuronal_sanger.obs['cell_type'] == leiden_from] = celltype
In [12]:
neuronal_harvard.obs['cell_states'].cat.categories
Out[12]:
In [14]:
trans_from=[['NRXN1_N1','NRXN4_N4'],['CD9+N0'],['Str2','Str5','Str6'],['doublets_N3'],['N7']]
trans_to = ['NRXN_Neuronal', 'CD9_Neuronal', 'Stromal-like', 'doublets', 'basic']
neuronal_harvard.obs['cell_type'] = [str(i) for i in neuronal_harvard.obs['cell_states']]
for leiden,celltype in zip(trans_from, trans_to):
for leiden_from in leiden:
neuronal_harvard.obs['cell_type'][neuronal_harvard.obs['cell_type'] == leiden_from] = celltype
Merge object¶
In [65]:
heart_neuronal = neuronal_harvard.concatenate(neuronal_sanger, neuronal_cells, neuronal_cd45, batch_key = 'cell_source', batch_categories = ['Harvard-Nuclei', 'Sanger-Nuclei', 'Sanger-Cells', 'Sanger-CD45'], join = 'outer')
heart_neuronal.shape
Out[65]:
In [66]:
pd.crosstab(heart_neuronal.obs['cell_source'], heart_neuronal.obs['cell_type'])
Out[66]:
In [67]:
heart_neuronal.obs['cell_type'] = heart_neuronal.obs['cell_type'].astype('category')
heart_neuronal.obs['cell_type'].cat.categories
Out[67]:
Remove immune-like stromal cells¶
In [68]:
heart_neuronal = heart_neuronal[~heart_neuronal.obs['cell_type'].isin(['doublets'])]
heart_neuronal.shape
Out[68]:
Normalise and scale dataset¶
In [69]:
neuronal_raw = heart_neuronal.copy()
sc.pp.normalize_per_cell(heart_neuronal)
sc.pp.log1p(heart_neuronal)
heart_neuronal.raw = heart_neuronal.copy()
In [70]:
sc.pp.highly_variable_genes(heart_neuronal, flavor = 'seurat', n_top_genes = 10000)
neuronal_hvg = heart_neuronal[:, heart_neuronal.var['highly_variable']]
sc.pp.scale(neuronal_hvg, max_value = 10)
Run scGEN on source¶
- Prepare the object
In [71]:
hvg_scGen_batch_1 = heart_neuronal.copy()
hvg_scGen_batch_1.obs["batch"] = heart_neuronal.obs["cell_source"].tolist()
hvg_scGen_batch_1.obs["cell_type"] = heart_neuronal.obs["cell_type"].tolist()
- Prepare scGEN object
In [72]:
tf.reset_default_graph()
network_batch_1 = scgen.VAEArith(x_dimension = hvg_scGen_batch_1.shape[1], model_path = "/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/")
- Train the model
In [73]:
network_batch_1.train(train_data = hvg_scGen_batch_1, n_epochs = 10)
- Correct batches
In [74]:
corrected_batch_1 = batch_removal_ct5(network_batch_1, hvg_scGen_batch_1)
- Visualise
In [75]:
sc.pp.neighbors(corrected_batch_1, random_state = 1786)
sc.tl.umap(corrected_batch_1, min_dist = 0.4, spread = 1, random_state = 1786)
sc.pl.umap(corrected_batch_1, color = ['cell_source','cell_type','region','donor','n_counts','n_genes','scrublet_score','gender'], size = 4, legend_fontsize = 6, color_map = 'RdPu', frameon = False)
Correct for donor¶
In [76]:
hvg_scGen_batch_2 = corrected_batch_1.copy()
hvg_scGen_batch_2.obs["batch"] = corrected_batch_1.obs["donor"].tolist()
hvg_scGen_batch_2.obs["cell_type"] = corrected_batch_1.obs["cell_type"].tolist()
- Create scGen model
In [77]:
tf.reset_default_graph()
network_batch_2 = scgen.VAEArith(x_dimension = hvg_scGen_batch_2.shape[1], model_path = "/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/")
In [78]:
network_batch_2.train(train_data = hvg_scGen_batch_2, n_epochs = 10)
In [79]:
corrected_batch_2 = batch_removal_ct5(network_batch_2, hvg_scGen_batch_2)
In [80]:
sc.pp.neighbors(corrected_batch_2, random_state = 1786)
sc.tl.umap(corrected_batch_2, min_dist = 0.4, spread = 2, random_state = 1786)
sc.pl.umap(corrected_batch_2, color = ['donor', 'cell_type', 'gender', 'cell_source', 'region', 'n_genes', 'n_counts', 'scrublet_score'], size = 6, legend_fontsize = 6, color_map = 'RdPu', frameon = False)
Call markers on Leiden¶
In [84]:
sc.tl.leiden(corrected_batch_2, resolution = 0.3, random_state = 1712)
sc.pl.umap(corrected_batch_2, color = ['leiden', 'cell_type'], size = 5, legend_fontsize = 6, color_map = 'RdPu', use_raw = True, frameon = False)
Call markers¶
In [85]:
sc.tl.rank_genes_groups(corrected_batch_2, 'leiden', method = 'wilcoxon', n_genes = 500, use_raw = True)
result = corrected_batch_2.uns['rank_genes_groups']
groups = result['names'].dtype.names
markers_immune = pd.DataFrame(
{group + '_' + key[:1]: result[key][group]
for group in groups for key in ['names', 'pvals_adj', 'logfoldchanges']})
markers_immune.head(10)
Out[85]:
Save annotated objects¶
In [86]:
heart_neuronal_raw = anndata.AnnData(X = neuronal_raw.X , obs = corrected_batch_2.obs, var = neuronal_raw.var, obsm = corrected_batch_2.obsm)
heart_neuronal_raw.raw = neuronal_raw.copy()
heart_neuronal_raw.write('/home/jovyan/experiments/2-subpopulation_manifolds/neuronal/2-aligned_batches/heart_neuronal_scGen_global_ctl200217.h5ad')
- Website